Gamma Function Subroutine

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;
}
Last edited on
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
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:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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.
Last edited on
Topic archived. No new replies allowed.