### C++ Subroutine Gamma Function

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
closed account (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.
 ``1234567891011121314`` ``````// 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:
 ``123456789101112131415161718`` ``````#include #include //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.
Last edited on
You could also test with the built-in gamma function from <cmath>

 ``12345678910`` ``````#include #include 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
closed account (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.
thank you for the help, this is actually much simpler than what I was trying to do.
closed account (D80DSL3A)
You're welcome. I got some useful review out of it too.
Topic archived. No new replies allowed.