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


implement Lanczos7;
 
include "sys.m"; sys: Sys;
include "draw.m";
include "math.m"; math: Math;
	lgamma, exp, pow, sqrt: import math;
 
Lanczos7: module {
	init: fn(nil: ref Draw->Context, nil: list of string);
};
 
init(nil: ref Draw->Context, nil: list of string)
{
	sys = load Sys Sys->PATH;
	math = load Math Math->PATH;
	# We ignore some floating point exceptions:
	math->FPcontrol(0, Math->OVFL|Math->UNFL);
	ns : list of real = -0.5 :: 0.1 :: 0.5 :: 1.0 :: 1.5 :: 2.0 :: 3.0 :: 10.0 :: 140.0 :: 170.0 :: nil;
 
	sys->print("%5s %24s %24s\n", "x", "math->lgamma", "lanczos7");
	while(ns != nil) {
		x := hd ns;
		ns = tl ns;
		# math->lgamma returns a tuple.
		(i, r) := lgamma(x);
		g := real i * exp(r);
		sys->print("%5.1f %24.16g %24.16g\n", x, g, lanczos7(x));
	}
}
 
lanczos7(z: real): real
{
	t := z + 6.5;
	x := 0.99999999999980993 +
		676.5203681218851/z -
		1259.1392167224028/(z+1.0) +
		771.32342877765313/(z+2.0) -
		176.61502916214059/(z+3.0) +
		12.507343278686905/(z+4.0) -
		0.13857109526572012/(z+5.0) +
		9.9843695780195716e-6/(z+6.0) +
		1.5056327351493116e-7/(z+7.0);
	return sqrt(2.0) * sqrt(Math->Pi) * pow(t, z - 0.5) * exp(-t) * x;
}
 