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


package require math
package require math::calculus
 
# gamma(1) and gamma(1.5)
 
set f 1.0
set f2 [expr {sqrt(acos(-1.))/2.}]
 
for {set x 1.0} {$x <= 10.0} {set x [expr {$x + 0.5}]} {
 
    # method 1 - numerical integration, Romberg's method, special
    #            case for an improper integral
 
    set g1 [math::calculus::romberg  \
                [list apply {{x t} {expr {$t ** ($x-1) * exp(-$t)}}} $x] \
                0 1 -relerror 1e-8]
    set g2 [math::calculus::romberg_infinity \
                [list apply {{x t} {expr {$t ** ($x-1) * exp(-$t)}}} $x] \
                1 Inf -relerror 1e-8]
    set gamma [expr {[lindex $g1 0] + [lindex $g2 0]}]
 
    # method 2 - library function
 
    set libgamma [expr {exp([math::ln_Gamma $x])}]
 
    # method 3 - special forms for integer and half-integer arguments
 
    if {$x == entier($x)} {
        puts [format {%4.1f %13.6f %13.6f %13.6f} $x $gamma $libgamma $f]
        set f [expr $f * $x]
    } else {
        puts [format {%4.1f %13.6f %13.6f %13.6f} $x $gamma $libgamma $f2]
        set f2 [expr $f2 * $x]
    }
}