Gamma Function Subroutine

[ShoxPhysics]

Hello,
I am new to working in C++ and in order to become better acquainted with it I am working on a program that I have no idea how to start. I am trying to create a subroutine that will calculate the gamma function of a real number, and the way that was suggested was to use the Weierstrass form. I have put the code below but it is not working properly, any suggestions would be greatly appreciated.

// 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;
}

[willf]

Well, I'm not sure if this is the programming issue, but the gamma function is continuous and therefore not always an interval. Try making it a double.

I'm new to programming my self, but I am likely starting grad school in physics this fall so I'm glad I'm not completely useless on this site! haha

[D80DSL3A]

You broke it! The function I gave here
http://www.cplusplus.com/forum/beginner/94174/
is a function of 2 variables which returns a double (as noted by willf).

...the way that was suggested was to use the Weierstrass form.

I didn't suggest that, you requested it!
From the previous thread:

I am trying to create a subroutine that will calculate the gamma function of a real number using the product series expansion of gamma.

That would be the Weierstrass form.

At any rate, I have improved on that function. It should be a function of 1 variable of course. The 2nd variable was for controlling precision.

Building a pre-determined precision into the function is easy:

double Gamma( double z )
{
    double g = 0.5772156649;// Euler-Mascheroni constant
    double retVal = z*exp(g*z);
    double x = 1.0;
    unsigned n=0;

    do
    {
        ++n;
        x = (1 + z/n)*exp(-z/n);
        retVal *= x;
    } while( x < (1 - 1e-12) );// the precision is determined here.

    return 1.0/retVal;
}

It gives decent results. The check values are from the std::tgamma() function as suggested by Cubbi in the previous thread.

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

But it's slow. The higher z is the longer it takes.
A better version takes advantage of the identity Gamma(z+1) = z*Gamma(z)

double Gamma_2( double z )
{
    double g = 0.5772156649;// Euler-Mascheroni constant
    double zi = floor(z);// integer part of z
    double zf = z - zi;// fractional part of z
    double retVal = zf*exp(g*zf);
    double x = 1.0;
    unsigned n=0;

    do  // calculate Gamma for the fractional part of z
    {
        ++n;
        x = (1 + zf/n)*exp(-zf/n);
        retVal *= x;
    } while( x < (1 - 1e-12) );

    retVal = 1.0/retVal;// invert

    // build up the value using the identity.
    for(double y = 0.0; y < zi - 0.5; ++y)
        retVal *= y + zf;

    return retVal;
}

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

That's almost 6x faster. If the time for only Gamma_2(5.1) was compared it would probably be larger than that, since that's the case that took the big time in the first case.

As you can see there is a small problem at integer z. The reason for the nan results can be seen at line 6: 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 10-15) can change it, so we arrive at line 17 with retVal still = 0;
Line 17 gives us the nan result (division by 0).
I'll let you modify the function to cope with the case of integer z.