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


let e = exp 1.
let pi = 4. *. atan 1.
let sqrttwopi = sqrt (2. *. pi)
 
module Lanczos = struct
  (* Lanczos method *)
  (* Coefficients used by the GNU Scientific Library *)
  let g = 7.
  let c = [|0.99999999999980993; 676.5203681218851; -1259.1392167224028;
	    771.32342877765313; -176.61502916214059; 12.507343278686905;
	    -0.13857109526572012; 9.9843695780195716e-6; 1.5056327351493116e-7|]
 
  let rec ag z d =
    if d = 0 then c.(0) +. ag z 1
    else if d < 8 then c.(d) /. (z +. float d) +. ag z (succ d)
    else c.(d) /. (z +. float d)
 
  let gamma z =
    let z = z -. 1. in
    let p = z +. g +. 0.5 in
    sqrttwopi *. p ** (z +. 0.5) *. exp (-. p) *. ag z 0
end
 
module Stirling = struct
  (* Stirling method *)
  let gamma z =
    sqrttwopi /. sqrt z *. (z /. e) ** z
 
end
 
module Stirling2 = struct
  (* Extended Stirling method seen in Abramowitz and Stegun *)
  let d = [|1./.12.; 1./.288.; -139./.51840.; -571./.2488320.|]
 
  let rec corr z x n =
    if n < Array.length d - 1 then d.(n) /. x +. corr z (x *. z) (succ n)
    else d.(n) /. x
 
  let gamma z = Stirling.gamma z *. (1. +. corr z z 0)
end
 
let mirror gma z =
  if z > 0.5 then gma z
  else pi /. sin (pi *. z) /. gma (1. -. z)
 
let _ =
  Printf.printf "z\t\tLanczos\t\tStirling\tStirling2\n";
  for i = 1 to 20 do
    let z = float i /. 10. in
    Printf.printf "%-10.8g\t%10.8e\t%10.8e\t%10.8e\n"
    		  z 
		  (mirror Lanczos.gamma z)
		  (mirror Stirling.gamma z)
		  (mirror Stirling2.gamma z)
  done;
  for i = 1 to 7 do
    let z = 10. *. float i in
    Printf.printf "%-10.8g\t%10.8e\t%10.8e\t%10.8e\n"
    		  z
		  (Lanczos.gamma z)
		  (Stirling.gamma z)
		  (Stirling2.gamma z)
  done