RosettaCodeData/Task/Gamma-function/Octave/gamma-function.octave

23 lines
604 B
Plaintext

function g = lacz_gamma(a, cg=7)
p = [ 0.99999999999980993, 676.5203681218851, -1259.1392167224028, \
771.32342877765313, -176.61502916214059, 12.507343278686905, \
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 ];
x=a;
if ( x < 0.5 )
g = pi / ( sin(pi*x) * lacz_gamma(1.0-x) );
else
x = x - 1.0;
t = p(1);
for i=1:(cg+1)
t = t + p(i+1)/(x+double(i));
endfor
w = x + double(cg) + 0.5;
g = sqrt(2.0*pi) * w**(x+0.5) * exp(-w) * t;
endif
endfunction
for i = 1:10
printf("%f %f\n", gamma(i/3.0), lacz_gamma(i/3.0));
endfor