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;
}