Introduction

In mathematics, when an exclamation mark, called a factorial, is used after a number it represents the multiplication of all the numbers from $0$ up to that number. For example, $6!=1\times 2\times 3\times 4\times 5\times 6\times=720$, so we can define the factorial as $n!=1\times 2\times \cdots n$. Now, the factorial function is not defined between integers, which is where the Gamma function comes in. The Gamma function, as oppose to the factorial, is a smooth function, meaning that it is defined between the integers and even for complex number. The relation between the factorial and the Gamma function is $$\Gamma(n)=(n-1)!$$ The Gamma function is defined through the improper integral $$\Gamma(n)=\int_0^\infty x^{n-1} e^{-x} dx$$ Note: The Gamma function is not the only function that interpolates the factorial. (link)
When I was looking into the Gamma function the integral seemed familiar, it looked similar to the extension of the Gaussian quadrature, which I saw in a textbook somewhere, called Gauss-Laguerre quadrature which is a numerical method for approximating improper integrals. According to a book by Adramowitz and Stegun, $$\int_0^\infty e^{-x} f(x)dx=\sum_{i=1}^n w_i f(x_i)+R_n$$ where $w_i=\frac{x_i}{(n+1)^2[L_{n+1}(x_i)]^2}$. I decided to use this to approximate the Gamma function.
But, first, we have to calculate the laguerre polynomials and from that the $i$-root of the laguerre polynomials then the weights from that. This needs only to be done once, so we can pre-calculate the weights and hard code them into the code we will write.

The Math

Since, we are only using one function with the Gauss-Lagurre quadrature, the intergral representing the Gamma Function, we don't need to calculate the weights multiple times but only once. If we compare the integral from the quadrature and the Gamma Function we can see that here $f(x)=x^{z-1}$, where we can set $z$ as a constant. Since this is only a demonstration I will only use 10-terms as to keep a relatively high amount of accuracy but still have a manageble number of terms to calculate.
We calculate the roots of the Laguerre polynomials, but first we calculate the Laguerre polynomials themselves. We can calculate the polynomials using the closed form formula $L_n(x)=\sum_{k=0}^n \binom{n}{k} \frac{(-1)^k}{k!}x^k$ Where we set $n = 10$ we get the polynomial $$1-10x+\frac{45}{2}x^2-20x^3+\frac{35}{4}x^4-\frac{21}{10}x^5+\frac{7}{24}-\frac{x^7}{42}+ \\ \frac{x^9}{896}- \frac{x^9}{36288}+\frac{x^{10}}{3628800}$$
For which, we then solve for the roots of this polynomial, in other words, when we set the polynomials equal to $0$. There are multiple ways of solving this polynomial, but I used a program I wrote which uses the Jacobi method, (link), but that is a topic for another post. Anyway, here are the roots
\begin{align}
x_1&=0.137793 \\
x_2&=0.729455 \\
x_3&=1.80834 \\
x_4&=3.40143 \\
x_5&=5.5525 \\
x_6&=8.33015 \\
x_7&=11.8438 \\
x_8&=16.2793 \\
x_9&=21.9966 \\
x_{10}&=29.9207
\end{align}
We now need to pass these rootes through the weights formula $w_i=\frac{x_i}{(n+1)^2[L_{n+1}(x_i)]^2}$,
\begin{align}
w_1&=0.308463 \\
w_2&=0.401115 \\
w_3&=0.218076 \\
w_4&=0.0620889 \\
w_5&=0.00950138 \\
w_6&=0.000753015 \\
w_7&=0.0000282585 \\
w_8&=4.24908\times 10^{-7} \\
w_9&=1.83954\times 10^{-9} \\
w_{10}&=9.91181\times 10^{-11}
\end{align}
Now that we have our roots we can set up the summation
$$\sum_{i=1}^n w_i f(x_i)=\sum_{i=1}^{10} w_i x_i^{z-1}\approx \Gamma(z)$$

The Code

We are now ready to implement the approximation of the Gamma function. Let's go ahead and do that in C. Let's create a file called gamma.c and include the following libraries we will need. Also, if compiling with gcc be sure to link the math library -lm.

#include <stdio.h>
#include <math.h>

Add the values we calculated as global static double arrays

static double x[10] = { 0.137793,
                        0.729455,
                        1.80834,
                        3.40143,
                        5.5525,
                        8.33015,
                        11.8438,
                        16.2793,
                        21.9966,
                        29.9207};

static double w[10] = { 0.308463,
                        0.401115,
                        0.218076,
                        6.20889e-2,
                        9.50138e-3,
                        7.53015e-4,
                        2.82585e-5,
                        4.24908e-7,
                        1.83954e-9,
                        9.91181e-11};

A for loop in a function will be used to carry out the summation

double gamma(double input)
{
  double s = 0;
  for (int i = 0; i < 10; i++)
    s += w[i] * pow(x[i], input - 1);
  return s;
}

In the main function, simply call the gamma function. To check the accuracy, call the function tgamma which is a function to calculate gamma provided by the math.h library.

int main()
{
  double param, result;
  param = 1.5;
  result = gamma(param);
  printf("Created gamma(%f) = %f\n", param, result);
  
  result = tgamma(param);
  printf("math.h tgamma(%f) = %f\n", param, result);
}

Here is the output of the program we just created.

ubrwlf@ubrpc:~/gamma$ gcc -o gamma gamma.c -lm
ubrwlf@ubrpc:~/gamma$ ./gamma
Created gamma(1.500000) = 0.889516
math.h tgamma(1.500000) = 0.886227

(Source Code)