// calculating a gamma function using a subroutine #include<iostream> #include<cmath> using std::cout; using std::endl; double PI = 3.141592654; int Gamma( int z, int Nterms ) { int g = 0.577216; int retVal = z*exp(g*z); // apply products for(unsigned n=1; n<=Nterms; ++n) retVal *= (1 + z/n)*exp(z/n); retVal = 1.0/retVal;// invert return retVal; } int main() { int z,gamma=1; cout<<"Enter the number to calculate gamma: "; cin>> num; for(int i=1;i<=num;i++) { gamma=Gamma(i); } cout<<"The gamma of the given number is: "<<gamma<<endl; cin.get(); return 0; } 
...the way that was suggested was to use the Weierstrass form. 
I am trying to create a subroutine that will calculate the gamma function of a real number using the product series expansion of gamma. 


Gamma(1/2) = 1.77245 check = 1.77245 Gamma(0.999) = 1.00058 check = 1.00058 Gamma(1) = 0.999999 check = 1 Gamma(3) = 2 check = 2 Gamma(5.1) = 27.9317 check = 27.9318 Process returned 0 (0x0) execution time : 1.781 s 


Gamma_2(1/2) = 1.77245 check = 1.77245 Gamma_2(0.999) = 1.00058 check = 1.00058 Gamma_2(1) = nan check = 1 Gamma_2(3) = nan check = 2 Gamma_2(5.1) = 27.9318 check = 27.9318 Process returned 0 (0x0) execution time : 0.328 s 
double retVal = zf*exp(g*zf);
where zf = fractional part of z, which = 0 in these cases. With our initial value for retVal = 0 no amount of multiplication (lines 1015) can change it, so we arrive at line 17 with retVal still = 0;