C++ Subroutine Gamma Function

[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 using the product series expansion of gamma. Any help would be greatly appreciated

[D80DSL3A]

Here's a simplistic implementation. Of course Nterms is supposed to be infinite, but that's not possible here. Try the function for various values of Nterms to see how rapidly it converges on the correct value.

// implements Weirstrass's form (infinite product)
double Gamma( double z, unsigned Nterms )
{
    double g = 0.577216;// Euler-Mascheroni constant
    double 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;
}

In all cases below I found that Nterms = 1,000 is needed to obtain accuracy to just 3 or 4 significant figures. It appears to converge slowly.
Testing with:

#include<iostream>
#include<cmath>
//using namespace std;
using std::cout;// better practice
using std::endl;

double PI = 3.141592654;

int main(void)
{
    cout << "Gamma(1/2) = " << Gamma(0.5, 1000) << endl;
    cout << "sqrt(PI) = " << sqrt(PI) << endl << endl;// checking value above
    cout << "Gamma(1) = " << Gamma(1.0, 1000) << endl;// 1
    cout << "Gamma(3) = " << Gamma(3.0, 1000) << endl;// 2
    cout << "Gamma(5.1) = " << Gamma(5.1, 1000) << endl;// 24 +
    cout << endl;
    return 0;
}

Output:

Gamma(1/2) = 1.77223
sqrt(PI) = 1.77245

Gamma(1) = 0.9995
Gamma(3) = 1.99103
Gamma(5.1) = 27.5716

EDIT: Source of formula: Mathematical methods for Physicists, 4th ed. page 594.
There are several other forms which may converge faster.

[Cubbi]

You could also test with the built-in gamma function from <cmath>

#include <iostream>
#include <cmath>
int main()
{
    std::cout << "Gamma(1/2) = " << std::tgamma(0.5) << '\n'
              << "sqrt(PI) = " << std::sqrt(std::acos(-1)) << '\n'
              << "Gamma(1) = " << std::tgamma(1.0) << '\n'
              << "Gamma(3) = " << std::tgamma(3.0) << '\n'
              << "Gamma(5.1) = " << std::tgamma(5.1) << '\n';
}

demo: http://ideone.com/ov4UzW

[D80DSL3A]

Thank you Cubbi. I was unaware that cmath supplies a gamma function. I assumed it would be too specialized so I didn't bother checking.

Obviously that's a better way of testing the function.

[ShoxPhysics]

thank you for the help, this is actually much simpler than what I was trying to do.

[D80DSL3A]

You're welcome. I got some useful review out of it too.