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


 
Module PrepareLambdaFunctions {
      Const e = 2.7182818284590452@
      Exp= Lambda e (x) -> e^x
      gammaStirling=lambda e (x As decimal)->Sqrt(2.0 * pi / x) * ((x / e) ^ x)
      Rad2Deg =Lambda pidivby180=pi/180 (RadAngle)->RadAngle / pidivby180
      Dim p(9)
      p(0)=0.99999999999980993@, 676.5203681218851@,   -1259.1392167224028@,  771.32342877765313@
      p(4)=-176.61502916214059@,  12.507343278686905@,  -0.13857109526572012@,  0.0000099843695780195716@
      p(8)=0.00000015056327351493116@
      gammaLanczos =Lambda p(), Rad2Deg, Exp (x As decimal) -> {
            Def Decimal a, t
            If x < 0.5 Then =pi / (Sin(Rad2Deg(pi * x)) *Lambda(1-x)) : Exit
            x -= 1@
            a=p(0)
            t = x + 7.5@
            For i= 1@ To 8@ {
                  a += p(i) / (x + i)
            }
             = Sqrt(2.0 * pi) * (t ^ (x + 0.5)) * Exp(-t) * a  
      }
      Push gammaStirling, gammaLanczos
}
Call PrepareLambdaFunctions
Read gammaLanczos, gammaStirling
Font "Courier New"
Form 120, 40
document doc$="     χ                        Stirling                     Lanczos"+{
}
Print $(2,20),"x", "Stirling",@(55),"Lanczos", $(0)
Print
For d = 0.1 To 2 step 0.1
      Print   $("0.00"), d,
      Print  $("0.000000000000000"), gammaStirling(d),
      Print  $("0.0000000000000000000000000000"), gammaLanczos(d)
      doc$=format$("{0:-10}  {1:-30}   {2:-34}",str$(d,"0.00"), str$(gammaStirling(d),"0.000000000000000"), str$(gammaLanczos(d),"0.0000000000000000000000000000"))+{
      }
Next d
Print $("")
clipboard doc$
 