46 lines
1.3 KiB
Plaintext
46 lines
1.3 KiB
Plaintext
/* From Rosetta Fortran */
|
|
test: procedure options (main);
|
|
|
|
declare i fixed binary;
|
|
|
|
on underflow ;
|
|
|
|
put skip list ('Lanczos', 'Builtin' );
|
|
do i = 1 to 10;
|
|
put skip list (lanczos_gamma(i/3.0q0), gamma(i/3.0q0) );
|
|
end;
|
|
|
|
|
|
lanczos_gamma: procedure (a) returns (float (18)) recursive;
|
|
declare a float (18);
|
|
declare pi float (18) value (3.14159265358979324E0);
|
|
declare cg fixed binary initial ( 7 );
|
|
|
|
/* these precomputed values are taken by the sample code in Wikipedia, */
|
|
/* and the sample itself takes them from the GNU Scientific Library */
|
|
declare p(0:8) float (18) static initial
|
|
( 0.99999999999980993e0, 676.5203681218851e0, -1259.1392167224028e0,
|
|
771.32342877765313e0, -176.61502916214059e0, 12.507343278686905e0,
|
|
-0.13857109526572012e0, 9.9843695780195716e-6, 1.5056327351493116e-7 );
|
|
|
|
declare ( t, w, x ) float (18);
|
|
declare i fixed binary;
|
|
|
|
x = a;
|
|
|
|
if x < 0.5 then
|
|
return ( pi / ( sin(pi*x) * lanczos_gamma(1.0-x) ) );
|
|
else
|
|
do;
|
|
x = x - 1.0;
|
|
t = p(0);
|
|
do i = 1 to cg+2;
|
|
t = t + p(i)/(x+i);
|
|
end;
|
|
w = x + float(cg) + 0.5;
|
|
return ( sqrt(2*pi) * w**(x+0.5) * exp(-w) * t );
|
|
end;
|
|
end lanczos_gamma;
|
|
|
|
end test;
|