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


fpprec: 30$
 
gamma_coeff(n) := block([a: makelist(1, n)],
   a[2]: bfloat(%gamma),
   for k from 3 thru n do
      a[k]: bfloat((sum((-1)^j * zeta(j) * a[k - j], j, 2, k - 1) - a[2] * a[k - 1]) / (1 - k * a[1])),
   a)$
 
poleval(a, x) := block([y: 0],
   for k from length(a) thru 1 step -1 do
      y: y * x + a[k],
   y)$
 
gc: gamma_coeff(20)$
 
gamma_approx(x) := block([y: 1],
   while x > 2 do (x: x - 1, y: y * x),
   y / (poleval(gc, x - 1)))$
 
gamma_approx(12.3b0) - gamma(12.3b0);
/* -9.25224705314470500985141176997b-15 */