problem with gamma computation

[Buck Evan]

The incomplete gamma function is resulting in negative outputs with positive real input, which shouldn't happen.

sage: gamma(60, 30).numerical_approx()

-1.28306738270893e74

Attached is a notebook showing the same in graphical form.

If I feed in numbers of a higher precision, I get the correct result:

sage: R10 = RealField(200)

sage: gamma(R10(60), R10(30)).numerical_approx()

1.38682990237888e80

My final problem is that there doesn't seem to be any way to get higher precision numbers to feed all the way through the plot function:

sage: show(plot3d(log(gamma(k, k*x)), (x, R10(0), R10(5/2)), (k, R10(1), R10(101))))

Launched jmol viewer for Graphics3d Object

(you would see a 3d graph with a bite out of it where the precision issue exists)