23 lines
603 B
Plaintext
23 lines
603 B
Plaintext
function x=gammal(z) // Lanczos'
|
|
lz=[ 1.000000000190015 ..
|
|
76.18009172947146 -86.50532032941677 24.01409824083091 ..
|
|
-1.231739572450155 1.208650973866179e-3 -5.395239384953129e-6 ]
|
|
if z < 0.5 then
|
|
x=%pi/sin(%pi*z)-gammal(1-z)
|
|
else
|
|
z=z-1.0
|
|
b=z+5.5
|
|
a=lz(1)
|
|
for i=1:6
|
|
a=a+(lz(i+1)/(z+i))
|
|
end
|
|
x=exp((log(sqrt(2*%pi))+log(a)-b)+log(b)*(z+0.5))
|
|
end
|
|
endfunction
|
|
|
|
printf("%4s %-9s %-9s\n","x","gamma(x)","gammal(x)")
|
|
for i=1:30
|
|
x=i/10
|
|
printf("%4.1f %9f %9f\n",x,gamma(x),gammal(x))
|
|
end
|