#include <math.h>
#include <conio.h>
#include <string.h>
#include <stdlib.h>
#include <stdio.h>

void integer_method(void);
void integer_check(void);
void zerocheck(void);
void limitcheck(void);
void gamma_calc(void);


const double pi=atan2(0.,-1.);


double a,abefore,aelk,u;

double lower, upper, ustore, ucheck, uu;

long i,k,m,n,nk;

int iflag = 0;
int ic;
int ipol=1;

void main()
{
	printf("\n\n\n  gamma,  ver 1.0,  \n\n  By Tom Irvine   Email: tomirvine@aol.com \n");
	printf("\n This program calculates the gamma function of a number. ");
	printf("\n The number may be any rational number.");



	while(1)
	{
        ic=0;
		iflag=0;
		ipol=1;

	
		printf("\n\n Enter number \n");

		scanf("%lf", &u);

		uu=u;

		zerocheck();
		

		if( fabs(u-double(0.)) >= 1.0e-06 )
        {
			limitcheck();

			u-=1;
	
			ic=1;

			ucheck = u;

			integer_check();

			if(ic==1 && iflag==0)
			{
				nk=1;

				printf("\n\n calculating.");
			
				ustore = u+1;
				
				u=long(u);
				integer_method();
				
				lower = a;
				upper = lower*(u+1);


				abefore = 1.;

					
                gamma_calc();

			}
					

			if( iflag == 0 )
			{

				printf("\n\n gamma function of %12.7g = %14.8g ",uu,aelk);
				
			} 

		}
	
		printf("\n\n");
		printf("\n press 1 for another calculation ");
		printf("\n       2 to exit \n");

		scanf("%ld", &i);

		if( i == 2 ){ break;}
	}

 
	
}
void integer_method(void)
{

		a=1.;

		n=long(u);

		for(i=1; i<=n; i++)
		{
			a*=i;
		}
}
void integer_check(void)
{


	int jk;

    for(jk=0; jk<=200; jk++)
    {


		if( fabs(ucheck-double(jk)) < 1.0e-06)
		{
			u=long(jk);

			integer_method();


			aelk=a;

			ic=0;

			break;
		}
        if( uu < 0.)
		{
			if( fabs(   fabs(uu)-double(jk)   ) < 1.0e-06)
			{
				printf("\n\n gamma function of %12.7g = infinity ",uu);
				iflag=1;
			}
		}
	}

}
void zerocheck(void)
{
	if( fabs(uu-double(0.)) < 1.0e-06 )
	{
		printf("\n\n gamma function of %12.7g = infinity ",uu);
	}

}
void limitcheck(void)
{
			if(fabs(u) > 170 )
			{
				printf("\n\n Number must be  <= 170 \n\n Press any key to continue.\n");
				getch();
				iflag=1;
		
			}

}
void gamma_calc(void)
{
		for(m=1; m<=8000000; m++)
		{
						if(nk==30000)
						{
							printf(".");
							nk=0;
						}
						nk++;
				
						double nx = pow(m,ustore);

						a=1.;

						if(m==1)
						{
							for(k=1; k<=m; k++)
							{
								a*=double(k);
								a/=(ustore + double(k));
							}
							a/=ustore;
						}
						else
						{
							a=abefore*double(m)/(ustore+double(m) );
						}


						abefore = a;

//						printf("\n m=%ld a=%lf",m,a );

						aelk=a*pow(m,ustore);
		}
}