Numerical methods for calculating derivative of gamma function

Question

I am working a mathematical statistic problem. It requires me to find MLE for a Gamma distribution by using numerical methods (rather than Method of Moments). Because Gamma distribution has two parameters, I have to solve a system of two non-linear equations which involve gamma function and its derivatives. I know there are methods that approximate gamma functions, Spouge's approximation or Lanczos approximation. (This post has a good summary Algorithm to compute Gamma function.) But I don't know how to approximate the derivative of a gamma function, that is, $$ \Gamma'(z) = \int_0^\infty (z-1) t^{z-1} e^{-t} \ln t \; dt . $$ I guess one can use Riemann sum to approximate this value, since the tail when $t \to \infty$ must be bounded. But this may be inaccurate. Is there any better way to approximate this integral (derivative of gamma function). Thanks a lot!