45 lines
1.1 KiB
Plaintext
45 lines
1.1 KiB
Plaintext
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;
|
|
}
|