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