Task

Implement one algorithm (or more) to compute the Gamma (



Γ


{\displaystyle \Gamma }

) function (in the real field only).

If your language has the function as built-in or you know a library which has it, compare your implementation's results with the results of the built-in/library function.

The Gamma function can be defined as:

  



Γ
(
x
)
=


∫

0


∞



t

x
−
1



e

−
t


d
t



{\displaystyle \Gamma (x)=\displaystyle \int _{0}^{\infty }t^{x-1}e^{-t}dt}


This suggests a straightforward (but inefficient) way of computing the 



Γ


{\displaystyle \Gamma }

 through numerical integration.



Better suggested methods:

 Lanczos approximation
 Stirling's approximation


#include <math.h>
#include <numbers>
#include <stdio.h>
#include <vector>
 
// Calculate the coefficients used by Spouge's approximation (based on the C
// implemetation)
std::vector<double> CalculateCoefficients(int numCoeff)
{
    std::vector<double> c(numCoeff);
    double k1_factrl = 1.0;
    c[0] = sqrt(2.0 * std::numbers::pi);
    for(size_t k=1; k < numCoeff; k++)
    {
        c[k] = exp(numCoeff-k) * pow(numCoeff-k, k-0.5) / k1_factrl;
        k1_factrl *= -(double)k;
    }
    return c;
}
 
// The Spouge approximation
double Gamma(const std::vector<double>& coeffs, double x)
{
        const size_t numCoeff = coeffs.size();
        double accm = coeffs[0];
        for(size_t k=1; k < numCoeff; k++)
        {
            accm += coeffs[k] / ( x + k );
        }
        accm *= exp(-(x+numCoeff)) * pow(x+numCoeff, x+0.5);
        return accm/x;
}
 
int main()
{
    // estimate the gamma function with 1, 4, and 10 coefficients
    const auto coeff1 = CalculateCoefficients(1);
    const auto coeff4 = CalculateCoefficients(4);
    const auto coeff10 = CalculateCoefficients(10);
 
    const auto inputs = std::vector<double>{
        0.001, 0.01, 0.1, 0.5, 1.0,
        1.461632145, // minimum of the gamma function
        2, 2.5, 3, 4, 5, 6, 7, 8, 9, 10, 50, 100, 
        150 // causes overflow for this implemetation
        };
 
    printf("%16s%16s%16s%16s%16s\n", "gamma( x ) =", "Spouge 1", "Spouge 4", "Spouge 10", "built-in");
    for(auto x : inputs) 
    {
        printf("gamma(%7.3f) = %16.10g %16.10g %16.10g %16.10g\n", 
            x,
            Gamma(coeff1, x),
            Gamma(coeff4, x), 
            Gamma(coeff10, x), 
            std::tgamma(x)); // built-in gamma function
    }
}
 