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


 
(import (scheme base)
        (scheme inexact)
        (scheme write))
 
(define PI 3.14159265358979323846264338327950)  
(define e 2.7182818284590452353602875)
 
(define gamma-lanczos
  (let ((p '(676.5203681218851 -1259.1392167224028 771.32342877765313 
             -176.61502916214059 12.507343278686905 -0.13857109526572012
             9.9843695780195716e-6 1.5056327351493116e-7)))
    (lambda (x)
      (if (< x 0.5)
        (/ PI (* (sin (* PI x)) (gamma-lanczos (- 1 x))))
        (let* ((x2 (- x 1))
               (t (+ x2 7 0.5))
               (a (do ((ps p (cdr ps))
                       (idx 0 (+ 1 idx))
                       (res 0.99999999999980993 (+ res 
                                                   (/ (car ps)
                                                      (+ x2 idx 1)))))
                    ((null? ps) res))))
          (* (sqrt (* 2 PI)) (expt t (+ x2 0.5)) (exp (- t)) a))))))
 
(define (gamma-stirling x)
  (* (sqrt (* 2 (/ PI x))) (expt (/ x e) x)))
 
(define gamma-taylor
  (let ((a (reverse
             '(1.00000000000000000000  0.57721566490153286061 
               -0.65587807152025388108 -0.04200263503409523553  
               0.16653861138229148950 -0.04219773455554433675
               -0.00962197152787697356  0.00721894324666309954 
               -0.00116516759185906511 -0.00021524167411495097  
               0.00012805028238811619 -0.00002013485478078824
               -0.00000125049348214267  0.00000113302723198170 
               -0.00000020563384169776 0.00000000611609510448  
               0.00000000500200764447 -0.00000000118127457049
               0.00000000010434267117 0.00000000000778226344 
               -0.00000000000369680562 0.00000000000051003703 
               -0.00000000000002058326 -0.00000000000000534812
               0.00000000000000122678 -0.00000000000000011813  
               0.00000000000000000119 0.00000000000000000141 
               -0.00000000000000000023  0.00000000000000000002))))
    (lambda (x)
      (let ((y (- x 1)))
        (do ((as a (cdr as))
             (res 0 (+ (car as) (* res y))))
          ((null? as) (/ 1 res)))))))
 
(do ((i 0.1 (+ i 0.1)))
  ((> i 2.01) )
  (display (string-append "Gamma ("
                          (number->string i)
                          "): "
                          "\n --- Lanczos : "
                          (number->string (gamma-lanczos i))
                          "\n --- Stirling: "
                          (number->string (gamma-stirling i))
                          "\n --- Taylor  : "
                          (number->string (gamma-taylor i))
                          "\n")))
 