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


/* From Rosetta Fortran */
test: procedure options (main);
 
  declare i fixed binary;
 
  on underflow ;
 
  put skip list ('Lanczos', 'Builtin' );
  do i = 1 to 10;
     put skip list (lanczos_gamma(i/3.0q0), gamma(i/3.0q0) );
  end;
 
 
lanczos_gamma: procedure (a) returns (float (18)) recursive;
    declare a float (18);
    declare pi float (18) value (3.14159265358979324E0);
    declare cg fixed binary initial ( 7 );
 
    /* these precomputed values are taken by the sample code in Wikipedia, */
    /* and the sample itself takes them from the GNU Scientific Library */
    declare p(0:8) float (18) static initial
         ( 0.99999999999980993e0, 676.5203681218851e0, -1259.1392167224028e0,
         771.32342877765313e0, -176.61502916214059e0, 12.507343278686905e0,
         -0.13857109526572012e0, 9.9843695780195716e-6, 1.5056327351493116e-7 );
 
    declare ( t, w, x ) float (18);
    declare i fixed binary;
 
    x = a;
 
    if x < 0.5 then
       return ( pi / ( sin(pi*x) * lanczos_gamma(1.0-x) ) );
    else
       do;
          x = x - 1.0;
          t = p(0);
          do i = 1 to cg+2;
             t = t + p(i)/(x+i);
          end;
          w = x + float(cg) + 0.5;
          return ( sqrt(2*pi) * w**(x+0.5) * exp(-w) * t );
       end;
  end lanczos_gamma;
 
end test;