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


/*
Here is the upper incomplete Gamma function. Omitting or setting
the second parameter to 0 we get the (complete) Gamma function.
The code is based on: "Computation of Special Functions" Zhang and Jin,
John Wiley and Sons, 1996
*/
 
SetFormat FloatFast, 0.9e 
 
Loop 10 
   MsgBox % GAMMA(A_Index/3) "`n" GAMMA(A_Index*10) 
 
GAMMA(a,x=0) {  ; upper incomplete gamma: Integral(t**(a-1)*e**-t, t = x..inf) 
   If (a > 171 || x < 0) 
      Return 2.e308   ; overflow 
 
   xam := x > 0 ? -x+a*ln(x) : 0 
   If (xam > 700) 
      Return 2.e308   ; overflow 
 
   If (x > 1+a) {     ; no need for gamma(a) 
      t0 := 0, k := 60 
      Loop 60 
          t0 := (k-a)/(1+k/(x+t0)), --k 
      Return exp(xam) / (x+t0) 
   } 
 
   r := 1, ga := 1.0  ; compute ga = gamma(a) ... 
   If (a = round(a))  ; if integer: factorial 
      If (a > 0) 
         Loop % a-1 
            ga *= A_Index 
      Else            ; negative integer 
         ga := 1.7976931348623157e+308 ; Dmax 
   Else {             ; not integer 
      If (abs(a) > 1) { 
         z := abs(a) 
         m := floor(z) 
         Loop %m% 
             r *= (z-A_Index) 
         z -= m 
      } 
      Else 
         z := a 
 
      gr := (((((((((((((((((((((((       0.14e-14 
          *z - 0.54e-14)             *z - 0.206e-13)          *z + 0.51e-12) 
          *z - 0.36968e-11)          *z + 0.77823e-11)        *z + 0.1043427e-9) 
          *z - 0.11812746e-8)        *z + 0.50020075e-8)      *z + 0.6116095e-8) 
          *z - 0.2056338417e-6)      *z + 0.1133027232e-5)    *z - 0.12504934821e-5) 
          *z - 0.201348547807e-4)    *z + 0.1280502823882e-3) *z - 0.2152416741149e-3) 
          *z - 0.11651675918591e-2)  *z + 0.7218943246663e-2) *z - 0.9621971527877e-2) 
          *z - 0.421977345555443e-1) *z + 0.1665386113822915) *z - 0.420026350340952e-1) 
          *z - 0.6558780715202538)   *z + 0.5772156649015329) *z + 1 
 
      ga := 1.0/(gr*z) * r 
      If (a < -1) 
         ga := -3.1415926535897931/(a*ga*sin(3.1415926535897931*a)) 
   } 
 
   If (x = 0)         ; complete gamma requested 
      Return ga 
 
   s := 1/a           ; here x <= 1+a 
   r := s 
   Loop 60 { 
      r *= x/(a+A_Index) 
      s += r 
      If (abs(r/s) < 1.e-15) 
         break 
   } 
   Return ga - exp(xam)*s 
}
 
/*
The 10 results shown:
2.678938535e+000  1.354117939e+000  1.0               8.929795115e-001  9.027452930e-001 
3.628800000e+005  1.216451004e+017  8.841761994e+030  2.039788208e+046  6.082818640e+062 
 
1.000000000e+000  1.190639348e+000  1.504575489e+000  2.000000000e+000  2.778158479e+000 
1.386831185e+080  1.711224524e+098  8.946182131e+116  1.650795516e+136  9.332621544e+155
*/